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We introduce ADAM, the All-Data Asteroid Modelling algorithm. ADAM is simple and universal since it handles all disk-resolved 
data types (adaptive optics or other images, interferometry, and range-Doppler radar data) in a uniform manner via the 2D Fourier 
transform, enabling fast convergence in model optimization. The resolved data can be combined with disk-integrated data (photome¬ 
try). In the reconstruction process, the difference between each data type is only a few code lines defining the particular generalized 
projection from 3D onto a 2D image plane. Occultation timings can be included as sparse silhouettes, and thermal infrared data are 
efficiently handled with an approximate algorithm that is sufficient in practice due to the dominance of the high-contrast (boundary) 
pixels over the low-contrast (interior) ones. This is of particular importance to the raw ALMA data that can be directly handled by 
ADAM without having to construct the standard image. We study the reliability of the inversion by using the independent shape 
supports of function series and control-point surfaces. When other data are lacking, one can carry out fast nonconvex lightcurve-only 
inversion, but any shape models resulting from it should only be taken as illustrative global-scale ones. 
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■ 1. Introduction 

O 

i Groundbased and other remote-sensing data on asteroids are ob- 
c/D tained with a variety of instruments that essentially sample re¬ 
gions on the surface of the target in various ways. These share 
some common mathematical characteristics of generalized pro¬ 
jections (Kaasalainen & Lamberg 2006; Kaasalainen 2011; Vi- 
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[ikinkoski & Kaasalainen|2014| . The most abundant source of 
data for asteroid shape and spin reconstruction is disk-integrated 
photometry, because e ven datasets sparse in time are often suf - 


ficient for modelling ( Kaasalainen| 


2004 


Durech et al. 


2006). 


Lightcurve-inversion procedures (Kaasalai nen et al.||20Ql 


are 


available at, e.g., the Database of Asteroid Models from Inver¬ 
sion Techniques (DAMIT) sit^] Due to the inherently limited 
information content of the disk-integrated data, the correspond¬ 
ing models arejisually most reliably described in convex space 


(Durech & Kaasalainen (2003), and further discussed below). 


However, even partially disk-resolved data offer a realistic pos¬ 
sibility of more detailed modelling. Previously described ap¬ 


proaches for such reconstruction are the SHAPE software (Ostro 
|et al.|2002|) for radar and lightcurve data, and the KOALA pro 


|et al . 2002) l or radar and lig htcurve data, and the KOALA 

cedure ( [Kaasalainen & Viikinkoski|20f2 |Carry et al.|2012 

optical images, occultation timings, and lightcurves. 


1 . 1 . The whole is more than the sum of its parts 

The best way to reconstruct a model of an asteroid is to use all 
available data. To combine disk-resolved data (adaptive optics 
or other images, interferometry, and range-Doppler radar data) 
with disk-integrated data (photometric or infrared lightcurves) 

1 http://astro.troj a.mff.cuni.cz/proj ects/asteroids3D 


and occultation timings (“sparse silhouettes”), we need a general 
procedure for using any data sources in asteroid modelling. We 
call this ADAM: All-Data Asteroid Modelling. Concise accounts 
of the various data types and their modelling aspects are given in 


Durech et al. 


and 

companion to those reviews. 


Kaasalainen & La mberg (2006), Kaasalainen & Durech (2013), 


(201). This paper is intended as a technical 


We present here the ADAM algorithm in a high-level format 
that includes all the necessary methods and formulae, either writ¬ 
ten here or given by references to the literature. We discuss and 
collect here the essential techniques and aspects of a complete 
inversion procedure capable of handling all the major asteroid 
data sources and formats. The key point is that complementary 
data sources can facilitate a good reconstruction even when none 
of them is sufficient alone. 


The paper is organized as follows. In Sect. 2 we describe 
the various shape supports we use in the reconstruction; some 
with the emphasis on global features, some concentrating on lo¬ 
cal details. This is intimately connected with the reliability esti¬ 
mate of the result, since independent shape representations help 
to reveal which features are the most probable ones. Sect. 3 in¬ 
troduces the Fourier transform method necessary for a simple 
and universal handling of data sources of disk-resolved type. In 
Sect. 4, we present examples of such types (interferometry, radar, 
and optical images). The interferometric data from ALMA are 
of particular interest here. We also discuss the special case of 
one-dimensional projections (continuous-wave radar and certain 
types of interferometry). In Sect. 5 we sum up everything in the 
form of the ADAM algorithm, and conclude in Sect. 6. Some 
basic ADAM functions are fisted in an Appendix. 
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1.2. ADAM software package 

Using the methods and algorithm described here and in 
Kaasalai nen et al.|(|2QQl| ), [Kaasalainen] ( |2011| ), and |Kaasalainen| 


& Viikinkoski (20 12[), writing an ADAM program from scratch 


is quite straightforward (for example, convex lightcurve inver¬ 
sion is inherently more complex). We have uploaded free-to- 
use ADAM code files and functions written in Matlab and C 
to a toolbox at the DAMIT site. These can be used for writ¬ 
ing customized inversion software, and for browsing and under¬ 
standing the computational methods. The latter, too numerous 
to be discussed here in detail, include techniques such as the 
partial derivative chains for gradient-based optimization, ray¬ 
tracing procedures, projections and transforms, scattering and 
luminosity models, GPU acceleration, etc. (Note that we do not 
offer any user support: the files are presented as is.) 

ADAM is a considerably more general package than the 
KOALA procedure (Kaasalain en & Vii kinkoski |2012[ |Carry| 
et al.||2012| ) that is based on extractable image contours. The 


KOALA contour-fitting principle is necessary for including oc- 
cultation data, so a full ADAM procedure inherits this function 
from KOALA. For fitting any pixel images, we recommend the 
ADAM Fourier-transform functions rather than KOALA. 

We take asteroid reconstruction to mean here that the fol¬ 
lowing output parameters are derived form input data: 1) shape 
(surface) definition, 2) rotational state (period and spin axis di¬ 
rection; possibly also terms for YORP acceleration, precession, 
or a binary orbit), 3) scattering or other luminosity parameters 
(often fixed a priori), and 4) image offset (alignment) and possi¬ 
ble other auxiliary or normalization terms. Without loss of gen¬ 
erality, we do not discuss each item separately, but mostly take 
the shape parameters to represent all the free ones since the 
optimization principle is technically the same for all parame¬ 
ter types. The speed, convergence, and reliability of gradient- 
based optimization methods are here superior to global methods 
(such as genetic algorithms or Monte Carlo; see the discussion in 
|Kaasalainen et al.|2001| ). We emphasize that the spin parameters, 
especially the period, usually have numerous local minima, so a 
dense enough comb of initial values of these is a prerequisite for 
a good final reconstruction. 


2. Shape 

Given the diverse shapes of asteroids and the continuing progress 
in instrument technology, effective methods for shape represen¬ 
tation are required for a general reconstruction scheme from ob¬ 
servations. In inverse problems it is typically not clear a priori 
how well a given shape support will perform. In this section we 
present shape supports and corresponding regularization func¬ 
tions well suited for asteroid-like shapes. 

2.1. Shape supports 

An important part of shape modelling is the choice of shape 
representation. Assuming a typical asteroid surface is homeo- 
morphic to the unit sphere, we can consider each coordinate 
as a function on the sphere, and choosing a suitable basis for 
functions, expand coordinate functions using this basis. This is 
straightforward to generalize to multiple bodies such as binaries. 
Typical such bases are spherical harmonics, spherical wavelets, 
and spherical splines. Our experiments suggest that parametriza- 
tions which expand each coordinate function separately tend to 
produce suboptimal results since they ignore the geometric de¬ 
pendencies and constraints between coordinates when consider¬ 
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ing surfaces represented by non-tangled meshes. Thus we have 
found it useful to consider two well-regulated but conceptually 
different shape supports in practice: octantoids based on spheri¬ 
cal harmonics, and subdivision surfaces. 

2.1.1. Function series 

An octantoid is a surface given by p e R 3 that can be 
parametrized in the form 


r x{o , (p) = e ^o,(p) s j n q cos ^ 

p{ 6 , (fi) = l y{ 6 , ip) = e a(0,ip)+b(9jp) s j n Q s j n ^ (1) 

[ z( 0 , ip) = e a(O,<p)+c( 0 ,<p) cos 0 ? 

where a , h and c are conveniently expressed as linear com¬ 
binations of the (real) spherical harmonic functions Y™(6,ip), 
with coefficients <z/ m , bi m and c/ m , respectively. Note that (0, ip ), 
0 < 6 < 7r, 0 < ip < 2n, are coordinates on the unit sphere 
S 2 parametrizing the surface but not describing any physical 
directions such as polar coordinates. As usual, the Laplace se¬ 
ries for a , b , c are useful for keeping the number of unknowns; 
i.e., the coefficients of TJ”, small and the surface smooth. If 
b = c = 0, this representation is the usual starlike one with 
the radius exp (a), but we have found that even if the target is 
starlike, the octantoid form allows the capture of detail better, 
and b and c can be represented with considerably fewer terms 
than the main function a. The number of shape parameters is 
thus between the (/ max + l) 2 of the starlike case and 3(/ max + l) 2 , 
when / max is the largest degree of the function series. The draw¬ 
back of this representation is its globality: one might want less 
smoothness regularization in some regions than in others. When 
more local control is desired (e.g., a feature clearly visible in fly¬ 
by images or in radar), the representation 0 may be expanded 
with spherical splines or spherical wavelets to provide local de¬ 
tail without affecting the global shape. Depending on the desired 
level of resolution and the non-starlike irregularity of the sur¬ 
face, the number of free function series coefficients is typically 
between 50 and 300 from low- to mid-resolution. Function series 
are seldom useful for high resolution, where one may ultimately 
want to adjust each vertex separately by defining individual a *, 
bi , and q. 

2.1.2. Subdivision control points 

Subdivision surfaces offer local control more than global repre¬ 
sentations like function series. Beginning with an initial set of 
vertices and corresponding triangles, called a control mesh, the 
surface is iteratively refined by adding new vertices and com¬ 
puting new positions for old vertices. The vertex coordinates of 
the control mesh form the parameter set defining the surface. 
Each subdivision step smoothes out the surface in a higher level 
of resolution. Well-behaving subdivision schemes converge to a 
smooth limit surface. 

In this paper we use the Loop subdivision scheme ( |Loop| 
|1987| . Considering a vertex p with immediate neighbours 
Po ,..., Pn-i > the subdivision method first creates new vertices 
by splitting each edge: 

3 P + 3 pi + 3pi-\ + p M 

q t = ---, i = 0, 1, (2) 

where the indices should be interpreted as modulo n. After the 
vertex creation step, the position of the vertex p is refined: 

p' = (1 - n[3)p + P^pi. ( 3 ) 
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The multiplier f3 is usually chosen to be 
5 (3 + 2cos(2 7i/n)) 2 


p= 1 - 

n 


8 


64 


( 4 ) 


but other choices are also possible. The limit surface is continu¬ 
ous; C 2 at the ordinary vertices (i.e. vertices that have 6 neigh¬ 
bours) and C 1 at extraordinary vertices. The number of free con¬ 
trol points for model rendering is similar to or somewhat lower 
than the number of function series coefficients (for a comparable 
level of resolution). 

The main computational aspect with subdivision methods is 
that the number of facets increases exponentially with the num¬ 
ber of divisions. After n subdivision steps, each facet that has 
been divided has produced 4 n subfacets. An alternative scheme 
to Loop subdivision is the V3-subdivision (Kobbelt 2000). In¬ 
stead of splitting the edges, V3-scheme subdivides facets by in¬ 
serting a new vertex to the facet centroid and connecting it to 
the vertices of the facet (Fig. [T]). The main attraction of the V3- 
scheme compared to the Loop subdivision is the slower increase 
(3 n ) of facets, while performing similarly in the limit. 

In practice, it is usually a good idea to choose the initial con¬ 
trol mesh to be an ellipsoid or a scaled convex surface obtained 
from lightcurve inversion, with a suitable number of vertices for 
the mesh. The number of subdivision steps should be chosen 
carefully: while each subdivision increases resolution and sta¬ 
bility by spreading the influence of each parameter to a larger 
number facets, the computational burden grows exponentially. 
Instead of subdividing all the facets, better performance may be 
obtained with adaptive subdivision, where only facets benefiting 
from increased resolution are subdivided. However, how to do 
this automatically during optimization is not obvious. A heuristic 
inclusion of surface regions to be refined based on a ranking of 
the improvement of the fit is one possibility (cf. the^ 2 -sensitivity 
map of [Kaasalainen & ViikinkoskT||2012| ); visual inspection of 
the model fit and a graphical user interface can guide in this. 


2.2. Regularization functions 

In inverse problems, finding a feasible regularization method is 
typically the most delicate part of problem solving. Ideally, both 
the shape representation and regularization method should be 
chosen to complement each other. The shape support should be 
general enough to represent probable shapes, and the regular¬ 
ization should prevent unrealistic or degenerate shapes while, at 
the same time, reveal the features present in the data. For octan- 
toids, the choice is remarkably easy. Assuming the basic shape 
is geometrically starlike, it is intuitively obvious to penalize the 
deviation from starlikeness. To this effect, we define 

n = L Z(fo L +c U (5) 

l,m 

Every starlike surface has a representation for which rj - 0, so 
rj is a natural quantity to be included in the final ^-function to 
be minimized (Sect. 5). The^ 2 -sum contains both the goodness- 
of-fit measure and the regularizing functions that represent prior 
assumptions and expectations of the solution. 

Subdivision surfaces have somewhat different smoothness 
properties in this regard. It is well known that the Loop sub¬ 
division converges to a smooth surface, so each subdivision step 
will produce a smoother result. However, it is computationally 
expensive to take a large number of subdivision steps. Therefore 
it is advantageous to combine a few, usually two or three, subdi¬ 
vision steps with mesh-based regularization methods. 



Fig. 1: Original control mesh (left) with 18 vertices (54 co¬ 
ordinates) as the shape parameter set and 32 facets, after two 
V3-subdivision steps (middle), and after four subdivision steps 
(right). 


While not strictly necessary, it is convenient to assume that 
the triangular mesh representing the shape forms a manifold. 
This assumption makes the checking of shadowing and illumi¬ 
nation both conceptually and computationally simpler. Thus it is 
imperative to avoid self-intersections, as they introduce errors to 
the fitting process. One approach is to explicitly check for in¬ 
tersecting facets and retriangulate if required. However, triangu¬ 
lation and intersection tests are costly, and usually optimization 
steps leading to self-intersections are suboptimal. A better ap¬ 
proach is to prevent self-intersections in the first place. 

Regularization based on dihedral angles penalizes large an¬ 
gles between adjacent facet normals; i.e., the regularization 
prefers planar regions. We thus want to minimize 

n = L Wi / 1 _ v i ' V A (6) 

tjer 


where T are the facets of the mesh, and is the unit normal vec¬ 
tor corresponding to the facet k. The sum is over all those facets 
j that are adjacent to the facet i, and the weights wij are usu¬ 
ally chosen to be unity. As a special case, we may suppress only 
concave features, obtaining convex regularization (Kaasalainen 
|& Viikinkoski| 2012 | ): 


72 


2 /A/ 


1 ^Ajd-vrvj), 


( 7 ) 


where A t is the area of the facet i and the sum is over those facets 
j that are adjacent to the facet i and tilted above its plane. 

To prevent degenerate facets and maintain a homogeneous 
mesh, it is advantageous to inhibit large variations in facet areas: 


73 =Yj,Ai-Af, (8) 

i 

where A is the mean facet area of the polyhedron. 

In practice, the regularization functions rj and y^ are suffi¬ 
cient for octantoid surfaces, while y 2 and 73 are useful for the 
subdivision surfaces. Unrealistically sharp angles can be pre¬ 
vented with y \, but too large a weight will inhibit convergence. 

In addition to geometric considerations, one can use regu¬ 
larization based on physical constraints, such as the requirement 
for the rotation axis to be close to the largest principal axis of 
the inertia tensor dKaasalainen[2011[|Kaasalainen & Viikinkoskif 
|20l2) . 

2.3. Reliability estimates 

The octantoid representation or the subdivision mesh tend to pro¬ 
duce aesthetically pleasing, “asteroid-like” surfaces, but it is not 
initially obvious which surface features of the model are actually 
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Fig. 2: Model of asteroid (41) Daphne from adaptive optics im- pig. 3: Asteroid (6489) Golevka reconstructed from disk- 
ages, reconstructed as a subdivision surface (left) and an octan- integrated photometry. From left to right: Convex, octantoid and 
toid (right). subdivision surface. 


present in the data, and which are the side effects of the shape 
support and the regularization used. Conventionally, Markov 
chain Monte Carlo (MCMC) methods are used to obtain a reli¬ 
ability estimate for the model parameters. However, in our case, 
modelling and s ystematic errors usually dominate ( Kaasal amen| 
& Durech|2006 ), rendering the MCMC approach inefficient and 


inaccurate because the error distribution is not known (it is cer¬ 
tainly not random Gaussian as in standard MCMC). 


Moreover, the posterior distribution of shape parameters 
from MCMC will not really tell anything about the reliability 
of the model with respect to data, but only about the distribution 
of the estimate within the adapted shape support. We have found 
that this results in an overoptimistic conception of the reliability 
of the result, simply because the acceptable shape results can¬ 
not be probed widely enough using one shape support only - the 
Monte Carlo procedure focuses on too small regions of shape 
variation for both computational and geometric reasons. In addi¬ 
tion, the computation of the model fit is time consuming if the 
data set and parameter space are large, making MCMC estima¬ 
tion computationally expensive. 


To circumvent these obstacles, we have found the follow¬ 
ing approach fast and robust in practice. Any real feature of the 
model based on the data should also be present if another, inde¬ 
pendent model type such as shape support is used. When model 
errors dominate, it is thus better to sample the “model space” 
within some x 1 than the ^ 2 -space with some fixed model. As 
an example, shape models of the asteroid Daphne from adaptive 
optics images and photometry ( [Viikinkoski & Kaasalainen|2014| ) 
using both the octantoid representation and subdivision surfaces 
are shown in Fig. [2] The models are quite similar and fit the data 
equally well, and their difference gives an idea of the real level 
of resolution. MCMC probing with either shape support leads to 
unrealistically small differences (insignificant compared to those 
in Fig. [2]). Even the shape-support test is likely to produce too 
optimistic reliability limits; the model error can be further en¬ 
larged by, e.g., introducing random fluctuations in the scattering 
properties over the surface. This principle could be developed 
into a meta-level Monte Carlo procedure that probes the space 
of possible model types using latent parameters. 


We conclude that shape sampling based on a fixed model 
type, no matter how diligently done with Monte Carlo or other 
methods, leads to overoptimistic resolution with artificial de¬ 
tails. A typical example of this is the radar model of the asteroid 
Itokawa that portrayed imaginary detail at the resolution level 
expected from the data while not capturing even the large-scale 
features. There was nothing wrong with the model fit to the data 
as such: the inverse problem was simply nonunique (or very un¬ 
stable) due to the restricted observing geometries and instrumen¬ 
tal projection (Sect. 4.2), but the constrained shape support of the 
program did not reveal this (Ost ro et al.|20 05 ; Nola n et al.|2014| ). 
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2.4. Inversion with photometry only 


Since ADAM utilizes photometric data in addition to disk- 
resolved data, we note that ADAM can be used to reconstruct 
a model using photometric data only (simply by using only the 
photometric fit function from the toolbox). This is easy and fast 
to do (and the shape rendering is even faster than the standard 
convex inversion of lightcurves), but the result is inevitably un¬ 
reliable: it is well known that even sizable nonconvex shape fea¬ 
tures require high solar phase angles to show in disk-integrated 


data (Kaasalainen et al 


Kaasalainen & Durech 


2001 


Durech & Kaasalainen 


2003 


2006). This can be probed with the shape 


reliability approach above. 

As an example, we show reconstructed shapes of the asteroid 
Golevka in Fig.[3j based on the data in Kaasalaine n et al.| ( |2001| ). 
Both the subdivision method and the octantoid-based model dis¬ 
play additional detail not seen in the convex model. However, the 
details are not supported by the data: the convex model gives at 
least as good a fit as the nonconvex ones, as is almost always the 
case with lightcurves (so far the only case of a better nonconvex 
model fit to photometry is that of the asteroid Eger in Durech 


|et al.|2012) . Indeed, with Golevka and other ground truth cases 
(maps from space probe missions), even the lightcurve fit with 
the correct shape and the scattering model assumed in inversion 
is not better than that with the co nvex model ( Kaasalai nen et al.| 


2001; Kaasalainen & Durech 2006). This underlines the fact that, 
due to systematic errors, any best-^ 2 solution with photometry 
only is likely to miss the correct details. 

While the convex model yields the best global agreement 
with the radar-based Golevka model (see the comparison in 
Figs. 3 and 4 in [Kaasalainen et al.||2002| ), the nonconvex ones 
portray much of the general sharpness and ruggedness of the 
body even though their details are not correct. The convex shape 
presents something of a softened error envelope within which 
numerous local shape variations are possible (as if the target 
were seen unfocused), while the nonconvex representations are 
samplings of those variations. Their details coincide neither with 
each other nor with those in the radar-based model, but they are 
useful as illustrations and for probing the potential shape options 
(cf. the nonconvex examples in |Kaasalainen et al.|2004| ). 


3. Fourier transform and information content 

As discussed in |Viikinkoski & Kaasalainen| ( [2014| ), the Fourier 
transform (FT) facilitates a natural interpretation for the pixel 
size as the maximum frequency present in the data, and makes 
it easy to incorporate the impulse response function of the 
imaging system. It also makes the optimization procedure fast 
and straightforward, without the cumbersome aspects related to 
pixellated image fields and binned model image distributions. 
The principle of the ADAM approach is to compare, instead of 
the images themselves, a set of FT samples (typically some thou- 
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sands depending on the level of resolution) from the model im¬ 
age with those of the data image, and to iterate until the best fit 
is found. This is described in Sect. 5. 

Letting T be the set of facets forming a model polyhedron 
and P a projection operator, the two-dimensional Fourier trans¬ 
form of a projected polyhedron in the (£, 77)-plane is 

T(u, v) = Yj ff B ' «£ >1) e~ 2m(u TT d//, (9) 

where B t is the luminosity value of the facet i, and the function 
/(£, 77 ) is unity if the point projected on (£, 77 ) is visible and zero 
otherwise. As shown in [ Viikinkoski & Kaasalainen| ( |2014| ), we 
obtain by Green’s theorem (dividing a facet into subfacets if nec¬ 
essary so that we may assume I is constant within each subfacet) 


T(u,v) = Y J B i '£ j I ij (u,v), (10) 

i j 

where 

1 (b — d)u - (a - c)v 
u ’ 4 7r 2 (u 2 + v 2 ) (a - c)u + (b - d)v 

x ^-Imiau+bv) _ e ~2m{cu+dv)^ q 

for the 7 -th boundary line segment (oriented counterclockwise) 
of the facet i, with the end points ( a , b) and (c, d). 

The summation over the interior edges of a projected poly¬ 
hedron can be reordered by noting that each polygon edge in 
the interior is shared by two polygons, so a new factor B can 
be taken to be the difference between the two B t , and the edge 
term is computed only once. This explicitly shows why most 
of the information in the image is indeed from the limb and 
shadow boundary curves discusse d in [Kaasalainen] ( |2011| ) and 
Kaasalai nen & Viikinkoski ( 2012] ). The values of B for interior 
triangle edges are usually close to zero (indeed, they vanish for 
the geometric scattering B t = const.), so most of the weight is on 
the boundary edges. In practice, this is confirmed by the similar 
results for, e.g., the asteroid Daphne obtained by KOALA and 
ADAM. There is little real information in the interior pixels of 
adaptive optics images, but on the other hand their errors do not 
distort the result either: the difference between the KOALA and 
ADAM models (for the same initial values and shape support) is 
negligible. 

The role of boundary information can be understood when 
comparing to the extreme case of lightcurve data: if we sum the 
pixel brightnesses over the image as in photometry, all the local 
shape information in the image is lost, so the remaining infor¬ 
mation is considerably more dependent on the light-scattering 
properties that are never very well known. But with images the 
boundary contrast is always the largest, so it is sufficient to 
have some kind of reasonable scattering (or thermal distribu¬ 
tion) model to account for the interior pixel contrasts. Indeed, the 
uniqueness theorems on the image, interferometry, occultation, 
or ra dar data are based on the robust boundary contour informa¬ 
tion (|Kaasalainen||2011[ [Kaasalainen & Viikinkoski||2012] |Vi-| 
ikinkoski & Kaasalain en|2014| ). With disk-integrated data only, 
Minkowski stability is luckily on our side when using convex 
models ( [Kaasalainen et al.|2001HKaasalainen et al.|2002] ). 


4. Data sources 

The versatility of the ADAM algorithm allows the handling of 
different data sources with only minor changes to the instrument- 
dependent part of the procedure (essentially just the definition of 




Fig. 4: The antenna locations of ALMA (left) and corresponding 
uv-plane visibilities (right). Images generated with the CASA 
software package. 


the instrumental projection plane and the adopted point-spread 
function). In this section, we present diverse examples of shape 
reconstruction with ADAM using both simulated and observed 
data. 


4. 1. Interferometry and ALMA 

The interferometric imaging method differs radically from a typ¬ 
ical telescope; instead of observing the sky brightness directly, 
the interferometer samples the Fourier transform of sky bright¬ 
ness. Each antenna pair of the interferometric array determines 
one sample on the Fourier plane. The maximum separation be¬ 
tween antennas determines the maximum attainable resolution. 
The interferometer most relevant to asteroid shape studies is the 
Atacama Large Millimeter Array (ALMA) in the Chilean desert. 
In its full configuration, the interferometer will be capable of 
observing at the resolution of a few milliarcseconds at the wave¬ 
length of 0.3 mm, corresponding to the separation of 16 km be¬ 
tween antennas. 

Given the brightness distribution /(£, 77 ) on the plane-of-sky, 
the visibility function is defined as the integral 


V(u, v) 



/(£ rj) e- 2m(u * +v,l) 


d^, 


( 12 ) 


which is a two-dimensional Fourier transform of the brightness 
distribution. Each antenna pair, corresponding to the projected 
baseline on the plane-of-sky, samples the visibility function. 
When the visibility function is sampled on a sufficiently dense 
set, the Fourier transform can be inverted to obtain the bright¬ 
ness distribution /(£, 77 ). But since the function V(u, v) is mea¬ 
sured only at a finite number of points, the observed visibility 
function is 


V(u, v) = F(u, v) V(u, v), 


(13) 


where F(u, v) is a sampling function corresponding to the sam¬ 
pled points on the ( u , v)-plane. Thus the obtained brightness dis¬ 
tribution is actually 

I&ri) = mTi)*I(g,T,), (14) 


i.e., a convolution of the true brightness distribution with the in¬ 
verse Fourier transform /(£, 77 ) of the sampling function. Deduc¬ 
ing the true brightness distribution I from the partially measured 
brightness / is an inverse problem and there are several iterative 
algorithms to infer I from /, see, e.g., Labeyrie et al. (2006). 

While the images obtained from the interferometer are in¬ 
formative, the great advantage with ADAM is that the algorithm 
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Fig. 5: Simulated, uncorrupted images with 5 mas pixel size (left 
column), observed dirty images generated with CASA (middle) 
and the reconstructed low-resolution shape model (right). Note 
that the middle-column images are not needed in inversion; we 
use the direct FT data instead. The images are what would be 
seen if the raw data were deconvolved for viewing purposes as 
is usually done for ALMA targets. The test shape model is from 
Ostro et al. ( 2000| ). 


& Kaasalainen|f20T4] ). The fast analytical computations are then 


efficient in the optimization. A simple thermophysical model is 
sufficient for shape reconstruction, as the most relevant informa¬ 
tion is contained in the boundary data, which are quite robust 
with respect to the thermal model used. This is in contrast to the 
disk-integrated thermal data that are more sensitive to both the 
surface properties and the thermal model. 

For thermal infrared imaging, ALMA facilitates asteroid ob¬ 
servations at resolution levels previously attained only by range- 
Doppler radar. To explore the possibilities of ALMA for shape 
modelling, we use the Common Astronomy Software Applica¬ 
tions (CASA) package developed by National Radio Astronom¬ 
ical Observatory (NRAO) to simulate observations. 

Consider a hypothetical asteroid with geocentric and helio¬ 
centric distances of 1.5 and 1 AU, respectively. The thermal flux 
is observed at the 350 GHz band, a frequency located in an at¬ 
mospheric window. There are 11 observation runs, each obser¬ 
vation lasting 50 s with ten-second integration time. We choose 
an antenna configuration providing approximate resolution of 10 
mas, a resolution which is well within the capabilities of ALMA. 
The antenna configuration and the corresponding uv-plane sam¬ 
pling pattern are shown in Fig. [4] Uncorrupted plane-of-sky im¬ 
ages, with a resolution of five milliarcseconds, are displayed in 
the column on the left in Fig. [5] We use the CASA software to 
add realistic atmospheric noise to the observations. The result¬ 
ing dirty images, which are obtained by assuming that the un¬ 
sampled frequencies are zero, are shown in the middle column. 
These images are provided for illustration purposes only, since 
ADAM uses the wv-plane samples directly. 

To test the ADAM reconstruction method, we use a low- 
resolution octantoid representation with 75 shape parameters. 
We also fit a scaling term, common to all observations. Usu¬ 
ally it is a good idea to use scaling specific to each observation, 
but in this case we know that all the simulated observations are 
done in similar conditions, so the common scaling term is jus¬ 
tified. The reconstructed shape is displayed in the right column 
in Fig. [5] with the same observation geometries as for the model 
images. The small-scale detail is lost, which is to be expected 
due to the added atmospheric noise and coarse instrument reso¬ 
lution. However, the bifurcated shape is well recovered despite 
the noisy data (note that we used ALMA data only). The com¬ 
putation time for this reconstruction was a few minutes. For real 
observations, complementary data are often provided by other 
observation methods, e.g., disk-integrated photometric data are 
almost always available. 


works directly with the values of the visibility function obtained 
from the instrument. This approach has several distinct advan¬ 
tages: 

- Sparse data may be used (e.g. interferometry with a few base¬ 
lines) 

- The distribution of antennas does not cause bias, since the 
Fourier transform is not inverted 

- Possible artefacts caused by the inversion process are 
avoided 

- The dependence between different observations is taken au¬ 
tomatically into account 


To obtain the luminosity values for the model surface (i.e., 
the brightness factor B t for each facet) in the infrared regime 
of ALMA, we can use the Fourier-series approximation of 
|Nesvorny & Vokrouhlicky[ ( |2008| ) as discussed in Viikinkoski 


4.2. Radar data 


The mathematical principles of the feasibility and uniqueness 
of the inversion of range-Doppler images are discussed in Vi- 
ikinkoski & Kaasalainen| ( |2014| ). Here we consider some prac¬ 
tical issues related to shape reconstruction. While other imag¬ 
ing methods rely on detecting the radiation of the sun that is 
reflected or re-radiated from the asteroid, radar provides its own 
illumination, making it possible to observe an asteroid regardless 
of the position of the sun. Moreover, in contrast to the visible 
or infrared wavelengths, the frequencies used by the radar are 
not significantly distorted by the atmosphere. Additionally, the 
properties of the waveform may be carefully controlled to reveal 
structural details on the surface of the asteroid. These properties 
make it possible to obtain data resolution down to ten meters or 
less for near-Earth asteroids, but this does not immediately trans¬ 
late to the same model resolution because of the inverse problem 
(cf. the Itokawa example in Sect. 2.3). 
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Range-Doppler radar resolves an object both in the range and 
in the line-of-sight velocity that translates to the Doppler shift of 
the reflected pulse. The frequency spectrum may be extracted 
by taking the fast Fourier transform of the pulses corresponding 
to a particular range gate. The actual hardware implementation 
and the signal processing are complicated as the detected signals 
are below the noise level of the instrument ( [Ostro et~aL||2002| ). 
Fortunately, the technical specifics are not required for the ac¬ 
tual shape reconstruction, since the radar performance may be 
modelled by the point-spread function of the system. 

The point p = ( x,y,z ) on the asteroid’s surface can be trans¬ 
ferred to the range-Doppler frame (r, D) by the linear mapping 

r = (xcos0 + ysin0)sin0 + zcos0 (15) 

D = co sin 6 (x sin 0 - y cos 0), (16) 

where co is the rotation rate of the asteroid around the z-axis 
and (6 , cp) are the spherical radar direction coordinates as seen 
from the asteroid. In this mapping, the range-Doppler radar im¬ 
age brightness L may be written as an integral over the asteroid 
surface S : 

Or, D) = JJ h r [r - /Cp)] h D [D - D(p )] B(p) I(p) dS, (17) 

where ho and h r are the point-spread functions of the radar sys¬ 
tem, corresponding to the Doppler-shifted frequency D and the 
range r, respectively. Here I is the visibility function, which 
is unity if the point is visible to the radar and zero otherwise. 
This form is similarly defined for all generalized projections 
( Kaasalainen & Lamberg||2006| . The mapping p —> (r, D) is 
unique, but its inverse is many-to-one, so the inherent informa¬ 
tion content of a range-Doppler image is considerably smaller 
than that of an optical image of similar resolution. Thus, while 
the nominal resolution provided by radar may be unmatched by 
any other instrument, the drawback of radar imaging is the diffi¬ 
culty of the interpretation of the images. 

The radar scattering function is given by B , which is usually 
a simple cosine law 

B(p) = C[p(p)]\ (18) 

where p is the cosine of the angle between the surface normal 
and the radar direction. The constants C and n measure the sur¬ 
face reflectivity and the specularity of scattering, respectively. 
The validity of Eq. ( [18] ) for modelling the microwave scattering 
from the asteroid’s surface is a rather convoluted question. While 
the cosine law is quite simplified, it should be noted that as the 
reflected wave is formed in a complicated manner by the surface 
material whose properties and roughness are usually unknown, 
fully realistic modelling of the reflected wave is not computa¬ 
tionally feasible. However, as in the other disk-resolved cases, 
most of the information is contained in the boundary contours 
and is thus independent of the scattering model used. 

Assuming the asteroid is modelled as a polyhedron with tri¬ 
angular facets T, the integral ( [17] ) may calculated separately for 
each facet, after projecting each triangle T t as a triangle PT t on 
the range-Doppler plane: 

L(r, D) = y BJi f h r (r - /) h D {D - D') d/ d D\ (19) 

t^t J pT i 


where we have assumed that the visibility I and the scattering 
law B are constant within a triangle. 



Fig. 6 : Mid-resolution shape model of the asteroid 2000 ET 70 
reconstructed from radar images. Viewing directions are from 
the positive x, y, and z axes, respectively. 


Taking the Fourier transform on both sides, applying the con¬ 
volution theorem, and writing Ti(u , v) for the sum over the edges 
of a Fourier transformed triangle as in Sect. 3, we obtain 


£(u,v) = Yj Bi li H r (u) H D (v) 7~i(u, v), 


( 20 ) 


T t eT 


where H r {u) and H D (v) are, respectively, the Fourier transforms 
of h r and ho- 

Like any images, radar plots are seldom correctly aligned 
in some reference frame due to the errors in the center of mass 
prediction, so the actual position of a radar image with respect to 
the two-dimensional projection of the model must be determined 
during the optimization. The task of image alignment is further 
complicated by the peculiar asymmetric structure of radar im¬ 
ages, especially the bright leading edge, other possible ridges 
of strong reflectivity, and the fading farthest-range pixels. If the 
alignment information is unknown, it is usually a good idea to 
fit the image offsets first to a fixed shape, obtaining better initial 
positions that can be used in the shape optimization. 

To demonstrate the reconstruction method, we make a fast 
ADAM model of the near-Earth asteroid 2000 ET 70 . Our goal is 
to get a quick first look at an initial model (to be refined at will). 
The asteroid was observed during February 2012 at Arecibo 
and Goldstone observatories using 2380 and 8560 MHz range- 
Doppler radars (Naidu et aL] |2013] ). The images obtained from 
Arecibo have a resolution of 15 m in range and 0.075 Hz in fre¬ 
quency. Goldstone images have a somewhat lower resolution, 15 
to 75 m and 1 Hz, respectively. Our goal is to produce medium- 
scale detail in the reconstructed shape, so a typical model choice 
is an octantoid with / max ~ 10 and around 1500 facets. Our exam¬ 
ple is “first-result oriented” on purpose, so we assume no infor¬ 
mation about the instrument-specific distortions, or more impor¬ 
tantly, knowledge about the point-spread functions determined 
by the instrument and the processing routines of the radar signal. 
Thus the point-spread function used in the shape reconstruction 
is simply the two-dimensional delta function. 

For each data image, we fit, in addition to the shape parame¬ 
ters, the offset with respect to the model center of mass and the 
reflectivity term in Eq. (TSJ). The reconstructed middle-resolution 
shape is shown in Fig. 6 [and the model fit to the data in Fig. [7] 
The shape model fits the boundary contours of the radar images 
satisfactorily, but there are some differences in the interior de¬ 
tails. This is a consequence of the parametrization and facet size 
chosen for reconstruction. The interior could be reproduced in 
greater detail by choosing a different parametrization, for exam¬ 
ple locally adaptive subdivision surfaces, or by refining the po¬ 
sitions of individual vertices. The model dimensions, shape fea¬ 
tures, and spin parameters agree with those published by Naidu 


et al.| ( [2013] ) (the spin parameters are identical except for a 2 C 
difference in the pole latitude, well within error limits). 

The main point of the initial low to middle resolution is that 
the speed of ADAM is considerable, and a detailed knowledge 
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of the instrument or the surface scattering physics is not needed, 
so one obtains a first model very fast by just feeding in the im¬ 
ages. The middle-resolution radar-based reconstruction (using 
82 radar images) was computed in less than an hour on a stan¬ 
dard laptop, and GPU programming can reduce the computation 
time significantly. This makes possible a broad sampling of the 
parameter space or real-time experimenting with various mod¬ 
els. Once a lower-resolution model has been adopted as the final 
frame, it is straightforward to refine it further. However, this re¬ 
quires accurate information about the point-spread and scatter 
functions. 

4.3. Adaptive optics and other images 

Model reconstruction from adaptive optics images in the Fourier 
approach is extensively covered in |Viikinkoski & Kaasalainen| 
( |2014| ), along with an example of the reconstruction of the main 
belt asteroid Daphne from adaptive optics images and photom¬ 
etry (Fig. [2|. Other imaging data may be incorporated into the 
framework using a similar approach. For instance, flyby images 
are, from the viewpoint of the reconstruction algorithm, concep¬ 
tually identical to the AO images. This is one of the attractions 
of ADAM: at the bare minimum, the user does not need to know 
anything about the images except their projection matrix and 
epochs. 

We note that the photometric data were actually not even 
needed in reconstructing Daphne (except for a better period 
value than with AO images only). The shape results with or with¬ 
out photometry are similar. This shows that even sparse AO data 
are well sufficient for modelling asteroid spin states and shapes 
in detail. 


4.4. One-dimensional projection operators 

In the regime between disk-integrated and disk-resolved ob¬ 
servations there are one-dimensional operators that project the 
plane-of-sky onto a line. Typical examples are the continuous- 
wave (CW) Doppler spectra that measure the distribution of 
the reflected power in frequency only, and the fine guidance 
sensors (FGS) onboard the Hubble Space Telescope, measur¬ 
ing the brightness distribution along an instrument axis. One¬ 
dimensional projections are seldom sufficient for actual shape 
reconstruction, but they may contain useful information about 
the object’s size or indications about the bifurcated structure 
( |Kaasalainen & ViikinkoskT||2012| ), and combined with other 
sources they facilitate shape inversion. 

In both examples, the measurement can be written in the 
form 


S(x) 


■/ 


/(£, rj) P(x - £ cos y - p sin y) dp, 


( 21 ) 


where I{£,i 7 ) is the plane-of-sky brightness (optical or radar) 
distribution of an object, P is the point-spread function of the 
instrument, and the angle y corresponds to the rotation of the 
sensor in the image plane. In Kaasalainen & Viikinkoski ( |2012| ), 
the integral was evaluated using a Monte Carlo method: the pro¬ 
jected model was sprinkled with uniformly distributed sampling 
points, and the integral was approximated as a sum over the vis¬ 
ible and illuminated sampling points. Here we demonstrate how 
the Fourier transform method can be used to interpret the integral 
as a tomographic operator on the Fourier plane. 

Taking the Fourier transform on both sides and using the 
projection-slice theorem (a slice of a 2-D FT along a line through 
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the origin equals the 1-D FT of the projection of the original 2- 
D function onto a line in the same direction; see e.g. Bracewell 
|2QQ3| ), we get 

S(f) = I if cos y, -/ sin y) p(f), ( 22 ) 

where the calligraphic characters denote the Fourier-transformed 
functions. Now it is obvious that S(f) is a slice of a Fourier- 
transformed brightness distribution along a line through the ori¬ 
gin, multiplied with the Fourier transform of the point-spread 
function. Moreover, this means that the same algorithm may be 
used to fit both FGS and adaptive optics data, and similarly both 
CW Doppler data and the range-Doppler images. In other words, 
we extract a one-dimensional Fourier transform from the 2-D 
model FT, and compare this with the 1-D FT formed from the 
data in the same way as in the full 2-D case. 


5. ADAM algorithm 

The flowchart in Fig. [ 8 ] describes the workings of ADAM. More 

specifically, the algorithm may be divided in five distinct steps: 

1. For each data image D { and observation geometry £*, the 
two-dimensional Fourier transform ( p'D l (u, v) of D t is sam¬ 
pled at a set of points {(m^-, v^)}, j = 1... N t , on the spatial 
frequency plane. The size of the set is chosen to correspond 
to the level of resolution. For pixel images, the transform 
can be computed by Eq. (T0| ) when considering each pixel as 
a polygon, or by using fast Fourier transform functions for 
chosen grid points (but the time spent for TDfu, v) is irrele¬ 
vant as most of the computations are for the trial models). 

2. The shape support and resolution level (number of parame¬ 
ters) are chosen. The parameters are initialized such that the 
initial shape is a sphere approximately equal in size to the 
target. 

3. For each observation geometry <5,, the Fourier transform 
TMfu, v) of the corresponding projection image M t of the 
model is calculated as described in the previous sections, 
together with the partial derivatives of TMfu, v) with re¬ 
spect to all optimized parameters. Ray-tracing, scattering or 
luminosity models, and coordinate transforms for the image 

[[Kaasalainen et al. ( 2001 


plane are discussed in 
( |2011| ), and | Viikinkoski & Kaasalainen| ( |2014| ). 


Kaasalainen 


4. An objective function ^ 2 is formed, with the square norm of 
the complex-valued FT fit error: 


vy) - Si(u t j, Vy) Vy 


i j= 1 


+ X>, 2 = : * 2 . 


(23) 

where (of of is the offset between the data image D t and the 
model image M t , and, by the convolution theorem, Si is the 
Fourier transform of the point-spread function of the imaging 
system. The y* represent various regularization terms defined 
above. 

For brevity, we have written only one data mode in Eq. ( [23] ); 
any number of modes with their goodness-of-fit functions 
can be added to the sum. These functions for photometry 
and silhouettes (occultations) are given in Kaasalainen et al. 
(2001), Kaasalainen (2011), and Viikinkoski & Kaasalainen 


(2014). The determination of the weights of the data modes 
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Fig. 7: Examples of range-Doppler images of the asteroid 2000 ET 70 from Arecibo observatory (rows 1 and 3) and corresponding 
simulated images from the mid-resolution model (rows 2 and 4). The contrast scale of the model image is somewhat modified to 
reveal inner image features. 


(as Ai for the regularization functions) is discussed in 
|Kaasalainen| ( | 201 1\ and |Kaasalainen & ViikinkosE| ( | 2012 | ). 
Weights can be determined for any subsets of data (e.g., less 
reliable images) if necessary. 

In addition, the intensity level of each data and model image 
must be normalized. Often it is enough to divide both model 
Mi and data image D t by their respective mean intensities. 
Equivalently, writing 

X 1 : = E \\£>i(uij, Vij ) - Miiuij, Vij) || 2 + Ay 2 , (24) 

ij 

we have 

ij 

where the mean (•) is taken over {(m^-, v^-)}, j = 1 ... N t . 
However, sometimes it is better to allow the intensity level 


-A'(^z j , Vi j ) Mi(Ui; , Vij) 


<IIAII> 


<Wi\\) 


+ Ay 2 , 


(25) 


of each M t to be a free parameter and use^ 2 ; this is useful 
in the case where the mean intensity of D t is corrupted by 
excessive noise in the image background (this is typical for 
range-Doppler images). This causes the ^ 2 el -based solution 
to have a slightly wrong size to compensate for the “diluted” 
normalized intensity level inside the actual object region of 

A- 

5. The shape and spin parameters and the offsets (« o *, c?) as well 
as the possible intensity level factors C; minimizing x 2 are 
determined with a suitable method such as the Levenberg- 
Marquardt algorithm. If there are several hundreds of pa¬ 
rameters, as in the case of fitting all shape vertices directly 
(instead of using function series or control points) to pro¬ 
duce maximal resolution, the conjugate gradient method is 
efficient ( |Kaasalainen et al.|2001| ). 
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Fig. 8: ADAM optimization algorithm as a schematic for one 
image type. 


6. Conclusions and discussion 

ADAM can handle radar data, images, interferometry (also in 
the thermal infrared), photometry, and occultations separately 
or in combinations. The ADAM procedure consists of a num¬ 
ber of modules, and there are various options for each mod¬ 
ule, customized to the end-user (e.g., the adopted optimization 
method, regularization functions, shape support and mesh struc¬ 
ture, ray-tracing method, coordinate system, luminosity/scatter 
model, image formats, etc.). In this sense, ADAM is a toolbox 
and a set of building blocks rather than a ready-made program. 

The main idea behind ADAM is the efficient use of the 
Fourier transform in handling both images and one-dimensional 
projection data. Fourier analysis has long been used in, e.g., im¬ 
age compression because it captures the essential information 
conveniently in a hierarchy of resolution. In the same vein, the 
FT approach in ADAM is ideal for producing models of desired 
levels of resolution, especially in the low- to medium-resolution 
category. In this framework, the goodness-of-fit function be¬ 
tween the model and the data is easy to compute and use in 
optimization. What is more, its convergence properties are more 
robust than if the images were used directly (in fact, one does not 
necessarily even have to look at the images or know much about 
the instrument that produced them). An analogy of this paradox 
is the simple one-dimensional problem of realigning two phase- 
shifted copies of a dual-frequency signal. If one does this by 
minimizing the signal difference by optimizing the shift in the 
original amplitude space, there are multiple local minima, but in 
frequency space the offset is found immediately. 

Despite its automatic character, ADAM should not be used 
as a black box: asteroid reconstruction is a complicated inverse 
problem, and one should be familiar with its mathematical prin¬ 
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ciples to understand the limitations and information content of 
the data sources. 
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Appendix: Sample ADAM functions 

In order to make the structure of ADAM more concrete, we show 
an example of how the program is divided into subroutines. We 
consider only the part of the program which computes the heat 
flux density of an object and its Jacobian (i.e., the partial deriva¬ 
tives of each modelled flux data point w.r.t. the free parameters), 
as all the other modules are structurally similar. In the case of op¬ 
tical images, the flux is replaced by brightness from scattering; 
for radar, by the signal strength in the range-Doppler plane. 

The partial derivatives w.r.t. the shape parameters are ini¬ 
tially calculated with respect to the vertex coordinates, making 
the routines independent of the parametrization used. Only in 
the final phase is the Jacobian determined using the chain rule. 
The functions are complex-valued, since the fitting is done on 
the frequency plane. In the optimization, the data are divided 
into real and imaginary parts and fitted separately. Usually, the 
Levenberg-Marquardt or the conjugate gradient method is used 
to optimize the^ 2 -fit. 

ADAM uses three different coordinates systems: the 
asteroid-centric coordinate frame with coordinate axes fixed to 
the asteroid, the asteroid-centric inertial frame, and the camera 
frame, which is determined by the instrumental orientation ge¬ 
ometry. The plane-of-sky view of an asteroid is obtained by pro¬ 
jecting the asteroid in the camera frame to the xy-plane. 

The Jacobian and the vector consisting of simulated values 
corresponding to the observations are computed using the fol¬ 
lowing subroutines: 

- Generate_HF_Matrix calls the subroutine 
Calc_Heat_Flux for each observation, and then com¬ 
bines the Jacobian submatrices into a full Jacobian matrix 
corresponding to all the observations. 

- Calc_Heat_Flux Calculates the Fourier transform of the 
two-dimensional plane-of-sky flux density and its par¬ 
tial derivatives by calling the subroutines Calc_Temp, 
Rot_Matrix, Cam_Matrix, Calc_Vis and Calc_FT. After 
the flux density of each facet is determined, the routine trans¬ 
forms the triangular mesh to the camera frame and projects 
the visible part of the mesh onto the xy-plane by discarding 
the z-coordinate. Finally, each vertex is transformed to the 
frequency plane using the Calc_FT subroutine and the con¬ 
tributions of the Fourier-transformed facets are summed. 

- Calc_Temp determines the temperature of the facets cor¬ 
responding to the observation geometry, using the FFT 
method. The partial derivatives of the temperature with re¬ 
spect to the shape parameters are also calculated. This sub¬ 
routine also calls subroutines Rot_Matrix and Calc_Vis. 

- Rot_Matrix calculates the rotation matrix needed to trans¬ 
form the object to the inertial frame. The rotation matrix is 
determined by the spin vector and the observation time. 

- Cam_Matrix determines the matrix needed to transform the 
inertial frame to the camera frame. This depends on the in¬ 
strument location and orientation. The z-coordinate codes 
the relative distance from the instrument. 

- Calc_Vis determines the visibility of facets using ray¬ 
tracing. In contrast to the optical case, a facet can be visible 
to the observer even if it is not illuminated by the sun. 

- Calc_FT Calculates the Fourier transform of a triangle pro¬ 
jected onto the xy-plane together with the corresponding par¬ 
tial derivatives. 

The most important setup factors determining the computa¬ 
tion time of ADAM are the numbers of facets and data points 
(the number of images and their pixels). The computation time 


increases approximately linearly with both numbers. The cost of 
actual optimization steps increases superlinearly with the num¬ 
ber of free parameters, but with large datasets (such as the radar 
example above) most of the computation is spent on determin¬ 
ing function values and their partial derivatives with respect to 
the vertex coordinates. In such cases, the number of shape pa¬ 
rameters is not critical for the computational cost in the mid¬ 
resolution regime, so one is free to choose a number that best 
corresponds to the resolution level (and set the number of facets 
accordingly). When the dataset is small, the computation time is 
short in any case, and the model is likely to be a low-resolution 
one, so again the number of parameters is not an issue. The 
cost of visibility determination by ray-tracing is insignificant as 
the potential blocker facets can be precomputed (Kaasalainen & 
|Torppa|2001| ). 

The shape reconstruction from observations is an easily par- 
allelizable problem. There are two obvious levels of parallelism: 
each observation can be calculated independently; or, within 
each observation, the contribution of each facet may be deter¬ 
mined simultaneously. The best choice depends on the com¬ 
puter’s architecture. The observation-level parallelism may be 
easily exploited using the MATLAB parallel computing tool¬ 
box, or more effectively by using the OpenMP API in the C 
language. The reduction in execution time scales almost linearly 
with the number of CPU cores. This is the approach currently 
implemented in ADAM. 

While it is possible to implement facet-level parallelism on 
the CPU by dividing the facet computations between several 
CPU cores, a more natural approach is to use one thread per 
facet. This kind of implementation is inefficient on the CPU, 
since the thread-switching latency is high compared to the run¬ 
ning time of a thread. However, the ability of the GPU to run 
thousands of lightweight threads simultaneously combined with 
the virtually costless thread switching makes it possible to attain 
orders of magnitude faster computation than with CPU. GPU- 
accelerated computing will be implemented in ADAM using the 
Nvidia CUDA programming platform. 


Article number, page 11 of[TT| 






